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Within the framework of density functional theory, the inclusion of exact exchange and non-local 
van der Waals/dispersion (vdW) interactions is crucial for predicting a microscopic structure of 
ambient liquid water that quantitatively agrees with experiment. In this work, we have used the 
local structure index (LSI) order parameter to analyze the local structure in such highly accurate 
ab initio liquid water. At ambient conditions, the LSI probability distribution, P(/), was unimodal 
with most water molecules characterized by more disordered high-density-like local environments. 
With thermal excitations removed, the resultant bimodal P(7) in the inherent potential energy 
surface (IPES) exhibited a 3:1 ratio between high- and low-density-like molecules, with the latter 
forming small connected clusters amid the predominant population. By considering the spatial 
correlations and hydrogen bond network topologies among water molecules with the same LSI 
identities, we demonstrate that the signatures of the experimentally observed low- (LDA) and high- 
density (HDA) amorphous phases of ice are present in the IPES of ambient liquid water. Analysis of 
the LSI autocorrelation function uncovered a persistence time of ~ 4 ps—a finding consistent with 
the fact that natural thermal fluctuations are responsible for transitions between these distinct yet 
transient local aqueous environments in ambient liquid water. 


I. INTRODUCTION 


Water is one of the most abundant molecules on Earth, 
yet the molecular arrangements in the solid, liquid, or 
even gas phases of water are extraordinarily complex. In 
the microscopic structure of water, this complexity arises 
due to a subtle balance between directional hydrogen- 
bonded interactions and weaker van der Waals (vdW) 
interactions which are non-directional 1- '. In fact, this set 
of intricate intermolecular interactions gives rise to many 
of the non-trivial thermodynamic and kinetic behaviors 
observed in water, particularly in the liquid phase, which 
significantly differs from many other simple liquids in 
properties such as the isothermal compressibility, heat 
capacity, and thermal expansivity, which all seem to di¬ 
verge at low temperature 8-10 . It has been argued that 
these anomalous properties are related to the observa¬ 
tion of two metastable glassy ice forms, the low-density 
(LDA) and high-density amorphous (HDA) ices, which 
can be prepared experimentally and are separated by a 
first-order-like phase boundary 11-16 . LDA and HDA re¬ 
covered at low temperature (~ 80 K) and ambient pres¬ 
sure differ by their densities, i.e., plda=0.93 g/cm 3 and 
Phda=1-17 g/cm 317 . The molecular environment in LDA 
is more similar to that of crystalline ice (cubic or hexag¬ 
onal) while that of HDA is more disordered and simi¬ 
lar to liquid water at standard temperature and pres¬ 
sure (STP). Furthermore, there is some experimental ev¬ 
idence for the existence of two metastable liquid states 
above the corresponding two glass transition tempera¬ 
tures 18-20 , and these two liquid extensions of the amor¬ 
phous ices are respectively called the low-density (LDL) 
and high-density liquids (HDL). In this regard, it has 


been suggested that the coexistence line separating LDA 
and HDA would extend into the liquid regime terminat¬ 
ing at a critical point called the liquid-liquid (LL) crit¬ 
ical point (LLCP) 21 , and this scenario would provide a 
possible explanation for the anomalous response proper¬ 
ties of water. However, it has not been possible to ex¬ 
perimentally prove two-liquid coexistence to date, since 
crystallization is too fast in macroscopic samples in the 
deeply supercooled regime. Very recently, however, wa¬ 
ter droplets with diameters of « 10 fin i were found to 
remain in the liquid phase for approximately 1 ms at 1 
bar and 227/^ K, i.e., at a temperature near the hy¬ 
pothesized LLCP 22 . The structure factors of these liquid 
droplets were measured with ultra-fast x-ray techniques 
and found to be very close to that of LDA, as one would 
expect for LDL 22 . While direct experimental evidence 
of a liquid-liquid coexistence in deeply supercooled wa¬ 
ter has been so far elusive, computer simulations have 
shown coexistence with certain classical model poten¬ 
tials (force fields) 23,24 such as ST2 25 , although other force 
fields, such as the coarse-grained mW model, do not show 
such behavior 26 . At temperatures above the hypothe¬ 
sized LLCP, undercooled water transforms continuously 
from LDL to HDL with increasing pressures as shown by 
neutron scattering data 2 '. Short-lived local fluctuations 
associated to LDL or HDL environments have long been 
suggested to show up in experimental observations of the 
liquid at STP 28-30 . 

To date, classification of the structures characteriz¬ 
ing the LDL and HDL molecular environments has been 
accomplished using classical force fields in conjunction 
with several order parameters to differentiate the local 
structural environments in these two forms of liquid wa- 
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ter 9,10,2 !’ 23 . In this regard, recent studies have shown 
that the local structure index (LSI) 31,32 is a very insight¬ 
ful order parameter that can be employed to characterize 
the LDL and HDL molecular environments 23 . For in¬ 
stance, the distribution of the LSI order parameter was 
found to be bimodal in the corresponding inherent poten¬ 
tial energy surface (IPES) of both supercooled and am¬ 
bient liquid water 33-35 . However, no attempt has been 
made yet to analyze the local environment in liquid water 
based on structures (configurations) generated using the 
more sophisticated ab initio based potentials available 
within the framework of density functional theory (DFT). 
While the typical structural relaxation times in the su¬ 
percooled regime are too long to be studied at present 
using DFT-based ab initio molecular dynamics (AIMD) 
simulations, these techniques can be utilized to accu¬ 
rately analyze the diverse local structural environments 
in ambient liquid water, which has also been the subject 
of considerable debate in the recent literature 28,30,36-38 . 

The predictive power of DFT-based AIMD simula¬ 
tions depends crucially upon the accuracy of the underly¬ 
ing exchange-correlation (XC) functional utilized in the 
quantum mechanical treatment of the electronic degrees 
of freedom. In this regard, a number of studies have now 
demonstrated that the most commonly utilized XC func¬ 
tionals based on the generalized gradient approximation 
(GGA) provide an inadequate description of liquid water 
at ambient conditions—a description in which the pre¬ 
dicted liquid is significantly overstructured and exhibits 
excessively sluggish dynamics 1,39-72 . These limitations of 
GGA-XC functionals stem from the well-known and dele¬ 
terious effects of self-interaction error' 3 as well as the ne¬ 
glect of the non-local electron correlation effects respon¬ 
sible for van der Waals (vdW) or dispersion interactions. 
As such, a reduction in the self-interaction error via the 
incorporation of a fraction of exact (Hartree-Fock) ex¬ 
change and an explicit inclusion of vdW/dispersion in¬ 
teractions have individually and collectively yielded a 
more accurate description of the structural and dynami¬ 
cal properties of liquid water 1 ’ 56,61 ’ 62 . 

In this regard, we have recently shown that uti¬ 
lization of the PBEO functional 74,75 , which includes 
25% exact exchange, in conjunction with a fully self- 
consistent (SC) implementation of the density-dependent 
vdW/dispersion correction of Tkatchenko and SchefHer 76 
(TS-vdW), i.e., the PBEO+TS-vdW(SC) XC functional, 
yields an oxygen-oxygen structure factor, Soo{Q), and 
corresponding radial distribution function, goo( r )i that 
are in quantitative agreement with the best available ex¬ 
perimental data 1 ’". This level of agreement between ab 
initio simulations and experiment primarily originates 
from an increase in the relative population of water 
molecules in the interstitial region between the first and 
second coordination shells, a collective reorganization in 
the liquid phase which is facilitated by a weakening of 
the hydrogen bond strength by the use of the PBEO hy¬ 
brid XC functional, coupled with a relative stabilization 
of the resultant disordered liquid water configurations 


by the inclusion of non-local vdW/dispersion interac¬ 
tions. This increasingly more accurate description of the 
underlying hydrogen bond network in liquid water also 
yielded other correlation functions, such as the oxygen- 
hydrogen radial distribution function, <?oh(?~), and the 
higher-order oxygen-oxygen-oxygen triplet angular dis¬ 
tribution, Pooo(@) i which encodes the degree of local 
tetrahedrality, as well as electrostatic properties, such as 
the effective molecular dipole moment, that are in bet¬ 
ter agreement with experiment 1 . Therefore, the over¬ 
all agreement between experiment and the PBEO+TS- 
vdW(SC) description of the microscopic structure of am¬ 
bient liquid water is indeed a very promising starting 
point for accurately measuring and further exploring the 
fluctuations in the local structural environments of liquid 
water by means of the LSI order parameter. 

In this work, we have extensively analyzed this con¬ 
tinuous distribution of local molecular environments in 
liquid water at ambient conditions and at the correspond¬ 
ing IPES using a systematic hierarchy of DFT XC func¬ 
tionals. For all of the XC functionals considered herein, 
the probability distributions of the LSI order parameter, 
P(I), were found to be unimodal in shape and rapidly 
decaying from maxima located in the low-LSI range, in¬ 
dicating that a majority of water molecules at ambient 
conditions are situated in more disordered high-density¬ 
like local environments. In this regard, we found that the 
proportion of molecules having low- and high-density-like 
environments varies significantly with the choice of the 
underlying XC functional, a trend which essentially man¬ 
ifests as a temperature effect. By removing thermal exci¬ 
tations and obtaining the IPES corresponding to ambient 
liquid water generated at the PBEO+TS-vdW(SC) level 
of theory, we found that the resultant P(/) was bimodal 
in shape and exhibited a 3:1 ratio between high- and 
low-density-likc molecules, with the latter forming small 
connected clusters amid the predominant population of 
high-density-like molecules. By considering the spatial 
correlations and the underlying hydrogen bond network 
topologies among IPES water molecules characterized by 
the same LSI identity, we demonstrate that distinctive 
signatures of the experimentally observed polymorphism 
in the amorphous phases of ice are also present in the 
IPES of ambient liquid water. Further analysis of the 
temporal correlations via the LSI autocorrelation func¬ 
tion has also uncovered a persistence time of ~ 4 ps—a 
finding which is consistent with the fact that LSI fluc¬ 
tuations in ambient liquid water occur due to natural 
thermal fluctuations between these distinct yet transient 
local aqueous environments. 

The remainder of the paper is organized as follows. In 
Section II, we describe the computational details of the 
AIMD simulations performed herein and the definition of 
the LSI order parameter. Section III contains an in-depth 
discussion of the results of these simulations and a com¬ 
parative analysis with the currently available theoretical 
and experimental literature. Finally, brief conclusions 
are provided in Section IV. 
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II. COMPUTATIONAL METHODS 

A. Simulation Details 

In this work, we have systematically performed a 
series of Car-Parrinello AIMD simulations 78 of am¬ 
bient liquid water usiug a hierarchy of different XC 
functionals. The sequence of XC functionals em¬ 
ployed herein includes the standard senri-local GGA 
of Perdew, Burke, and Ernzerhof (PBE) 79 , the cor¬ 
responding hybrid PBEO' 4,75 which includes 25% ex¬ 
act exchange, and the self-consistent (SC) dispersion- 
corrected analogs 80 thereof, i.e., PBE+TS-vdW(SC) 
and PBEO+TS-vdW(SC), based on the Tkatchenko- 
Scheffler 76 density-dependent vdW/dispersion func¬ 
tional. 

All of these AIMD simulations were performed in the 
canonical ( NVT ) ensemble using periodic simple cubic 
simulation cells with lattice parameter set to reproduce 
the experimental density of liquid water at ambient con¬ 
ditions. All of the AIMD simulations were initially equi¬ 
librated for approximately 2 ps and then continued for at 
least an additional 20 ps for data collection. Four AIMD 
simulations were performed on (H 2 0)64 at 300 K using 
the PBE, PBEO, PBE+TS-vdW(SC), and PBEO+TS- 
vdW(SC) XC functionals. Since a classical treatment of 
the nuclear degrees of freedom is insufficient for a quanti¬ 
tatively accurate description of the microscopic structure 
of ambient liquid water, we have also performed an ad¬ 
ditional AIMD simulation at 330 K on (H 2 0)i28 at the 
PBE0+TS-vdW(SC) level, a technique suggested by the 
lowest-order perturbative expansion in h of the free en¬ 
ergy to account for the quantum mechanical nature of the 
nuclear degrees of freedom 81 . In practice, this increase 
of approximately 30 K in the simulation temperature has 
been found to mimic the nuclear quantum effects (NQE) 
in structural quantities such as the oxygen-oxygen radial 
distribution function ( goo{f )) in both DFT 82 and force 
field 83,84 based MD simulations. 

All calculations reported herein were performed within 
the plane-wave and pseudopotential framework and uti¬ 
lized a modified development version of the Quantum 
ESPRESSO (QE) software package 85 . To meet the ad¬ 
ditional computational demands associated with large- 
scale AIMD simulations based on hybrid XC functionals, 
we have employed a linear scaling O (TV) exact exchange 
algorithm that exploits the natural sparsity associated 
with the real-space maximally localized Wannier function 
(MLWF) 86 representation of the occupied Kohn-Sham 
electronic states, which has been developed 87 and exten¬ 
sively optimized 88 in our research group. In addition, we 
have also developed and utilized a linear scaling O (N) 
self-consistent implementation of the TS-vdW dispersion 
correction 80 , which provides a framework for computing 
atomic Cq dispersion coefficients as explicit functionals 
of the charge density, i.e., Cq^ab = C(i,An[p{ Y )\-> thereby 
accounting for the local chemical environment of each 
atom' 6 . More explicit descriptions of the simulation de¬ 


tails and theoretical methods employed herein are pro¬ 
vided in Ref. 1 . 

For comparative purposes, we have also performed sev¬ 
eral simulations of liquid water using the TIP4P/Ice 89 
rigid force field, a classical water potential chosen be¬ 
cause of its ability to reproduce the experimental freezing 
temperature of liquid water to within 1 K. Using the DL- 
POLY code 90 , all classical simulations were performed for 
2-10 ns (after equilibration for 50-100 ps) on (H 2 0)256 
in the NVT ensemble with periodic simple cubic sim¬ 
ulation cells set to reproduce the experimental density 
of ambient liquid water. To control the ionic tempera¬ 
ture, Nose-Hoover chain thermostats 91 were employed in 
conjunction with an integration time step of 1 fs. The 
Lennard-Jones (LJ) potential was truncated at 9.5 A and 
all electrostatic contributions were computed using the 
standard Ewald summation technique 92 . 


B. Local Structure Index (LSI) 

To quantify the degree of inhomogeneity in the lo¬ 
cal molecular environments of liquid water, we have uti¬ 
lized an order parameter introduced by Shiratani and Sa- 
sai 31,32 , which associates a local structure index (LSI) to 
the individual water molecules comprising the liquid. In 
essence, the LSI order parameter is the mean-squared- 
deviation among the radial distances corresponding to 
the set of molecules that surround a given water molecule. 
To compute the LSI value for a given water molecule, the 
set of radial oxygen-oxygen distances {rj} corresponding 
to the N neighboring molecules that are within a cutoff 
distance of 3.7 A from the reference molecule are first 
ordered as follows: rq < r 2 < • ■ • < rj < rj +i < ■ • ■ < 
rjv < 3.7 A < rjv+i- The LSI value (/) is then defined as 
the inhomogeneity in this distribution of radial distances, 
i.e., 

1 N 

/= ^E [a ^u-<a>] 2 , (1) 

3 = 1 

in which (A) is the mean over all Aj+ij = Vj+\ — rj. 
Hence, I provides a convenient quantitative measure of 
the fluctuations in the distance distribution surround¬ 
ing a given water molecule within a sphere defined by a 
radius of ~ 3.7 A. In doing so, I measures the extent 
to which a given water molecule is surrounded by well- 
defined first and second coordination shells. 

The schematic diagram provided in Figure 1 illustrates 
the different local molecular environments that can be 
distinguished and quantified by the LSI order parameter. 
For instance, a molecule characterized by a high I value 
will be found in a more ordered local environment, in 
which the neighboring water molecules are concentrated 
in the region between 2.7 3.2 A and sparse in the region 
between 3.2-3. 8 A. This leads to a more prominent sep¬ 
aration between the first and second coordination shells 
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FIG. 1. Schematic description of high-density-like and lo¬ 
cally disordered (left) vs. low-density-like and locally ordered 
(right) environments, which correspond to low and high val¬ 
ues of the LSI order parameter, I, respectively. The dark and 
light blue areas represent the first and second coordination 
shells around the central water molecule, respectively. The 
separation between the water molecules in the first and sec¬ 
ond coordination shells (depicted by red and black circles, re¬ 
spectively) are less prominent in aqueous environments char¬ 
acterized by a low LSI value. 


and a relatively low local atomic number density (as de¬ 
picted on the right of Figure 1). On the other hand, 
a molecule characterized by a low I value implies that 
the molecule is situated in a locally disordered environ¬ 
ment, a consequence of which is a relatively high packing 
of neighboring molecules in the interstitial region and 
hence an increase in the local atomic number density (as 
depicted on the left of Figure 1). 


III. RESULTS AND DISCUSSION 

A. LSI Distribution in Ambient Liquid Water 

The diversity of local microscopic environments found 
in ambient liquid water can be characterized by the prob¬ 
ability density distribution of the LSI order parameter 
(P(/)) and we begin our discussion by analyzing how 
P(I) depends on the underlying DFT XC potential and 
the simulation temperature. As shown in Figure 2(a), the 
P (I) in ambient liquid water obtained from various DFT- 
based AIMD simulations are all unimodal in shape and 
rapidly decaying from maxima located between 0.013 A 2 
(PBE0+TS-vdW(SC) at 330 K) and 0.023 A 2 (PBE at 
300 K). This finding is indicative of the fact that a ma¬ 
jority of the molecules are situated in high-density-like 
locally disordered environments in ambient liquid water. 
To illustrate that the smaller LSI values indeed corre¬ 
spond to more disordered local environments, we have 
computed P(/) for the tetrahedrally-ordered ice I h phase 
at 300 K, which is quite broadly distributed around the 
relatively larger LSI value of 0.25 A 2 and vanishes at I = 
0.05 A 2 (Figure 2(a)). Considering the most accurate 


XC functional employed herein, PBE0+TS-vdW(SC), 
approximately 80% of water molecules were found with 
/ < 0.05 A 2 at 330 K—a simulation temperature that 
has been elevated by 30 K to approximately account for 
NQE in ambient liquid water. In this regard, one can im¬ 
mediately notice that a 30 K decrease in the simulation 
temperature (while keeping the underlying XC functional 
fixed) leads to a decrease of approximately 8% in the 
number of water molecules that are situated in disordered 
high-density-like local environments (i.e., as character¬ 
ized by I < 0.05 A 2 ). This is accompanied by an increase 
in the relative population of water molecules with more 
ordered low-density-like local environments, i.e., water 
molecules characterized by larger LSI values. 

To further explore the temperature dependence of 
P (I), we have employed the classical TIP4P/Ice poten¬ 
tial, a rigid water force field which is able to reproduce 
the experimental freezing point of liquid water to within 1 
K 89 . As shown in Figure 2(b), by systematically lowering 
the simulation temperature from 330 K to 210 K in incre¬ 
ments of 30 K (fixing the simulation cell corresponding 
to the experimental density of ambient liquid water), we 
observed a systematic decrease in the relative population 
of water molecules with lower LSI values (i.e., molecules 
in locally disordered environments) coupled with an in¬ 
crease in the relative population of water molecules with 
higher LSI values (i.e., molecules in locally ordered envi¬ 
ronments). Similar to the ab initio PBE0+TS-vdW(SC) 
simulations discussed above, the relative population of 
water molecules with / < 0.05 A 2 decreases by approx¬ 
imately 5% when the simulation temperature of the liq¬ 
uid is reduced from 330 K to 300 K using the TIP4P/Ice 
force held. Here, we note in passing that the depen¬ 
dence of P(/) on temperature has also been demonstrated 
with other classical water force fields 33-35 , yielding simi¬ 
lar qualitative results to those observed in this work. 

Quite remarkably, the dependence of P(I) on the 
underlying DFT XC potential essentially manifests as 
a temperature effect. As shown in Figure 2(a), the 
LSI distributions obtained with the PBE and PBE+TS- 
vdW(SC) XC functionals are markedly different from 
that obtained with PBE0+TS-vdW(SC), despite the fact 
that the simulation temperature was set to 300 K in 
all of these cases. For example, the peak height in 
the P (I) computed with PBE is « 43% smaller than 
that obtained at the PBE0+TS-vdW(SC) level of the¬ 
ory at 300 K. Consequently, only 50% of the PBE wa¬ 
ter molecules were characterized by I < 0.05 A 2 as op¬ 
posed to 72% in liquid water generated with PBE0+TS- 
vdW(SC) at 300 K. These differences clearly indicate 
that the molecules in PBE-generated liquid water sam¬ 
ple high-density-like disordered local environments to a 
much lesser extent than in liquid water obtained using 
more accurate vdW-inclusive hybrid XC functionals. As 
such, this enhanced population of molecules having more 
ordered low-density-like local structures in ambient PBE 
liquid water can readily explain the overstructured pair 
correlation functions and exceedingly sluggish dynamics 
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FIG. 2. Probability density distributions of the local structure index (P(I)) in liquid water obtained from (a) DFT-based AIMD 
simulations and (b) classical force field (TIP4P/Ice) based simulations. 


routinely observed at the PBE level of theory 40-42,45,46,55 . 

To follow up on this point, it is noteworthy to reiterate 
the fact that the variations in P (/) with the choice of the 
underlying XC potential show a significant resemblance 
to the systematic reductions observed in the peak heights 
of P(/) while supercooling liquid water using the classical 
TIP4P/Ice potential. In this regard, ab initio simulations 
employing the more accurate PBEO+TS-vdW(SC) XC 
functional at 330 K show a remarkable similarity with 
the TIP4P/Ice simulation at w 330 K 93 , whereas the 
P (/) obtained at the PBE level of theory more closely re¬ 
sembles the P (/) distributions generated by the classical 
TIP4P/Ice simulations performed between 210-240 K. 
From this observation, one can infer that PBE liquid wa¬ 
ter simulated at ambient temperature (~ 300 K) behaves 
more like a deeply supercooled liquid, which is effectively 
60-90 K lower than the actual simulation temperature, 
a finding which is in line with the previous independent 
estimate of the melting temperature of ice with the PBE 
XC functional, found to be approximately 100 K higher 
than the experimental value 55 . 

Over the past decade, liquid water generated by the 
PBE XC functional at ambient conditions has been re¬ 
ferred to as “supercooled” due to the following character¬ 
istic features, which include overstructuring 42 , underes¬ 
timated diffusion coefficient 46 , and overestimated melt¬ 
ing temperature 55 . From the analysis of P(/) above, the 
sluggishness associated with PBE liquid water can physi¬ 
cally be explained by this relatively increased population 
of low-density-like molecules—a population which is en¬ 
ergetically stabilized by the favorable hydrogen bonding 
motifs that are facilitated by the presence of a locally 
ordered environment. Since the fluidity in liquid water 
arises from the fact that the hydrogen bonds compris¬ 
ing the underlying tetrahedral network are continuously 
breaking and forming, the relatively large number of in¬ 
tact hydrogen bonds in this population of low-density¬ 


like water molecules leads to an exceedingly sluggish liq¬ 
uid characterized by a significantly underestimated dif¬ 
fusion coefficient. On the other hand, the formation and 
stabilization of an increased population of high-density¬ 
like water molecules—as found in the AIMD simulations 
that employ vdW-inclusive hybrid XC functionals—are 
primarily driven by entropic effects 94 , which lead to a 
decrease in the average number of intact hydrogen bonds 
per water molecule and therefore a significantly more dif¬ 
fusive liquid. This increasing disparity between liquid 
water generated at the PBE and PBE0+TS-vdW(SC) 
levels of theory will become even more apparent when the 
LSI distribution is computed from the inherent potential 
energy surface (IPES) 95-9 ' of ambient liquid water, i.e., 
when the thermal fluctuations present in the liquid are 
removed, which brings us to the subject of the next sec¬ 
tion. 


B. LSI Distribution in the Inherent Potential 
Energy Surface (IPES) of Ambient Liquid Water 

A collection of local potential energy minima result¬ 
ing from systematic quenches along a liquid trajectory, 
the IPES is a powerful concept that allows us to sepa¬ 
rate packing from thermal motion effects 95,96 . To study 
the LSI distribution in the corresponding IPES of ambi¬ 
ent liquid water, approximately 100 configuration snap¬ 
shots were taken at 0.25 ps intervals from the PBE (300 
K) and PBE0+TS-vdW(SC) (330 K) simulations and 
rapidly quenched to the nearest local potential energy 
minima via second-order damped dynamics 98 . Upon re¬ 
moval of the thermal energy from these ambient liquid 
water structures, the LSI distributions were found to be 
bimodal in shape (as opposed to the unimodal distribu¬ 
tions found above in Section III A) and are depicted in 
Figure 3. For both of these XC functionals, the two peaks 
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FIG. 3. Probability density distributions of the local structure 
index (P(/)) in the inherent potential energy surfaces (IPES) 
of liquid water obtained from DFT-based AIMD simulations 
and classical force field (TIP4P/Ice) based simulations 
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in the bimodal P (!) distributions are separated by distin¬ 
guishable minima (Ani n ) located at approximately 0.14- 
0.16 A 2 (which should be contrasted against the fact that 
Anin = 0.13-0.14 A 2 for the TIP4P/Ice force field over 
the temperature range of 210-330 K). Considering the 
fact that our IPES P(/) have a resolution of > 0.01 A 2 
(due to the relatively low number of uncorrelated config¬ 
urations accessed in our AIMD simulation), the position 
of this so-called isobestic point seems to be slightly de¬ 
pendent on these rather distinct ab initio and classical 
potentials. In this regard, previous studies of the IPES in 
liquid water based on classical force fields have reported 
that the position of this minimum is nearly invariant over 
a wide range of temperatures and pressures 33 ^ 35 . How¬ 
ever, we note in passing that the position of the isobestic 
point is also weakly dependent on the choice of the cutoff 
distance used in Equation (1) to define the LSI distribu¬ 
tion of radial distances 34 . 

In the IPES, we found that the relative heights of the 
two peaks in the bimodal P(/) distribution obtained with 
the PBE XC functional are significantly different than 
those obtained at the PBE0+TS-vdW(SC) level of the¬ 
ory. Based on the positions of Anj n = 0.14 A 299 in the 
IPES P(/) distributions, the relative amount of highl¬ 
and low-density-like molecules in liquid water generated 
by the PBE0+TS-vdW(SC) XC functional was found as 
77% and 23%, respectively, whereas in the case of PBE 
they are almost equally distributed. Hence, with thermal 


excitations removed, PBE0+TS-vdW(SC) liquid water is 
still predominantly comprised of water molecules in high- 
density-like locally disordered environments, while PBE 
liquid water is essentially an equal mixture of low- and 
high-LSI sites. As observed above in Section III A for 
the thermally excited P(/) distributions, the IPES P(/) 
obtained from PBEO-l-TS-vdW(SC) closely resembles the 
TIP4P/Ice classical force field results obtained at ~ 300 
K, whereas the PBE distribution is in closer agreement to 
the TIP4P/Ice data obtained at ~ 210 K. This finding is 
again suggestive that the IPES of liquid water simulated 
at ambient temperature with the PBE XC functional is 
more similar to the IPES of a deeply supercooled liquid. 

This P (/) analysis at the IPES again exemplifies the 
importance of exact exchange and vdW/dispersion forces 
in generating an accurate microscopic description of am¬ 
bient liquid water. The increased population of high- 
density-like molecular environments found in PBE0+TS- 
vdW(SC) liquid water at ambient conditions and at the 
corresponding IPES originates from the collective effects 
of weakening the directional hydrogen bonds (facilitated 
by a reduction in the self-interaction error) and strength¬ 
ening the non-directional vdW/dispersion interactions. 
Since this vdW-inclusive hybrid XC functional is able 
to describe ambient liquid water with good accuracy 1 , 
we perform the next set of analyses at the PBE0+TS- 
vdW(SC) level of theory only , in order to further our 
understanding of the spatial and temporal correlations 
among molecules identified by low- and high-density-like 
environments in liquid water in the IPES of ambient liq¬ 
uid water. 


C. Signature of Polyamorphism in the IPES of 
Ambient Liquid Water 

The bimodality observed in the IPES of liquid wa¬ 
ter allows us to characterize molecules having high- and 
low-density-like local structural environments based on 
the position of the isobestic point, I m i n = 0.14 A 299 , 
by dividing water molecules according to I < Ani n and 
I > Anin, respectively (see Section IIIB). In the IPES of 
liquid water simulated with PBE0+TS-vdW(SC) at 330 
K, the ratio between high- and low-density-like molecules 
is found to be approximately 3:1 based on this distinc¬ 
tion. When the effects of thermal motions are removed 
in the IPES, low-density-like molecules form small con¬ 
nected clusters amid the larger relative population of 
high-density-like molecules (see Figure 4). In our sim¬ 
ulation cell, which contains 128 water molecules, we ob¬ 
tained 29±7 low-density-like molecules in the IPES on 
average and found cluster sizes among them that con¬ 
tained up to 34 molecules. In this case, a molecule is 
considered to be part of a cluster if it is located within 
3.33 A of any other water molecule present in the clus¬ 
ter (based on the oxygen-oxygen distance). The clus¬ 
ters formed among the low-density-like water molecules 
do not span the length of the simulation cell. However, 
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FIG. 4. Snapshots of typical liquid water structures (con¬ 
figurations) found in the IPES of ambient liquid water com¬ 
puted with the PBEO+TS-vdW(SC) XC potential at 330K. 
The blue spheres represent water molecules characterized by 
I < 7min (i.e., high-density-like water molecules), whereas the 
red spheres represent water molecules characterized by I > 
7min (i.e., low-density-like water molecules). The green box 
refers to the simulation cell containing 128 water molecules 
with only the oxygen atoms shown above for clarity. 


these small connected clusters may act as precursors for 
the growth of larger clusters of low-density-like molecules 
upon further cooling, a trend which was found in other 
theoretical studies 33-35 . 

In order to quantitatively study how molecules char¬ 
acterized by these distinct local environments are spa¬ 
tially correlated in the IPES, we have computed oxygen- 
oxygen radial distribution functions, goo(i"), separately 
among the high- and low-density-like water molecules. 
In doing so, we found a remarkable similarity between 
the goo(r ) computed among the high-density-like (i.e., 
with / < / m in) water molecules only and the low-density¬ 
like (i.e., with / > / m in) water molecules only with 
the experimentally observed goo(r) for the high-density 
(HDA) and low-density (LDA) amorphous ice phases, re¬ 
spectively. In Figure 5(a)-(b), these theoretical goo(r) 
are compared against the goo ( r ) corresponding to the 
low-temperature (80 K) HDA and LDA phases, which 
were obtained from empirical potential structure refine¬ 
ment (EPSR) of experimental neutron diffraction data 17 . 
From this picture, two different molecular environments 
are clearly evident - -in the case of I < / m ; n (see Fig¬ 
ure 5(a)), the interstitial region between the first and 
second coordination shells is largely populated, whereas 
the separation between the first and second coordina¬ 
tion shells is prominent in the case of I > I min (see Fig¬ 
ure 5(b)). Beyond the first coordination shell, i.e., in the 
intermediate range order (IRO) region, correspondence 
between the theoretical and experimental goo ( r ) is re- 
markable up to the length of the simulation cell employed 


herein. The prominent feature of a much shallower and 
broader second shell present in HDA is observed among 
IPES molecules with / < / m i n (see Figure 5(a)), whereas 
a much sharper and narrower second shell at ~ 4.5 A 
present in LDA is observed among IPES molecules with 
/ > / m i n (see Figure 5(b)). The correspondence between 
theoretical and experimental goo ( r ) in the region of the 
first coordination shell, which defines the short range or¬ 
der (SRO), is significantly effected by the thermal and 
quantal motions present in HDA and LDA. The com¬ 
plete removal of thermal effects in the quenched struc¬ 
tures obtained at the IPES is reflected in the fact that 
the first peaks of both theoretical goo( r ) are too sharp 
when compared to the experimental goo (f) of the amor¬ 
phous structures, which show a much broader first peak. 
Due to the very same reason, we also observed a decrease 
in the population of molecules in the interstitial regions of 
the theoretical goo( r ) when compared to the analogous 
experimental quantities. As such, these findings clearly 
suggest that the IPES of ambient liquid water contains 
the signatures of the experimentally observed polymor¬ 
phism in the amorphous phases of ice. 

Apart from the distinctions observed in the goo(r) 
above, the structural differences between HDA and LDA 
also lie in the topologies of their underlying hydrogen 
bond networks. These differential topologies can be 
quantified using ring statistics analysis—a tool which 
has been instrumental in theoretically characterizing 
and delineating the evolution of the HDA and LDA 
phases 100,101 . In the case of LDA, which is more ordered, 
one finds a narrower distribution of rings centered about 
the most predominant six-membered ring, whereas in the 
case of HDA, which is more disordered, one finds a much 
broader distribution of rings with a substantial popula¬ 
tion having more than six members. In fact, the existence 
of larger rings provides a pathway of having more com¬ 
pact and disordered packing of water molecules in the 
amorphous ice and liquid phases of water 100,101 . Here, 
we have used a geometric definition of intact hydrogen 
bonds 102 in conjunction with the criteria of King 103 to 
find the shortest-path closed rings containing N mem¬ 
bers among water molecules characterized by I > / m i n 
only and I < I min only, and computed the resultant nor¬ 
malized probability distributions, P(iV), depicted in Fig¬ 
ure 5(c). From this figure, it can readily be seen that 
the hydrogen bond networks connecting water molecules 
with the same LSI identities have very different topolo¬ 
gies. For instance, we found a broad distribution of rings 
and a tail that reaches out to N = 12 for the high-density- 
like molecules. Closed rings containing N > 12 members 
were omitted because these rings start to include peri¬ 
odic images of the molecules forming the ring due to the 
finite size of the simulation cells employed in this work. 
On the other hand, the P (N) obtained from low-density¬ 
like molecules is much narrower and does not contain 
rings with more than N = 9 members, consistent with 
the observation that low-density clusters do not span 
the length of the cell. As such, these difference in the 








N 


FIG. 5. Comparison of the oxygen-oxygen radial distribution 
functions, goo{r), computed among water molecules charac¬ 
terized by (a) I < /min only and (b) I > /min only based 
on liquid water structures (configurations) taken from the 
IPES generated with the PBEO+TS-vdW(SC) XC potential, 
with the experimental structures of the low-temperature (80 
K) high-density (HDA) and low-density (LDA) amorphous 
ice phases, respectively. The experimental goo(r) were ob¬ 
tained using empirical potential structure refinement (EPSR) 
of neutron diffraction data 17 . Probability distributions of the 
hydrogen-bonded IV-membered rings, P(iV), computed from 
structures on the PBE0+TS-vdW(SC) IPES are depicted in 
panel (c). All P(N) were normalized to unity and therefore 
do not reflect the total number of rings of a given size. 


P (N) of low- and high-density-likc clusters bears strong 
similarities with those found for the LDA and HDA ice 
phases 100 ’ 101 . 

Instead of a predominance of N = 6 membered rings, 
as expected in the structures of LDA 100,101 , we found 
that the largest majority of rings contain N = 5 mem¬ 
bers with significantly large populations of smaller rings 
among the high- and low-density-like water molecules in 
the IPES of ambient liquid water. As such, this higher 
population of shorter rings (containing N < 6 members) 
is reflective of the fact that the IPES of ambient liq¬ 
uid water does not contain true phases of either HDA 
or LDA. In order to obtain the HDA and LDA phases 
of ice, structural relaxations on much larger length and 
time scales are required—relaxations which are certainly 
not achieved with instantaneous quenching to the near¬ 
est local potential energy minimum, the technique em¬ 
ployed in this work to generate the corresponding IPES. 
However, the topological differences in the underlying hy¬ 
drogen bond networks observed here is striking and con¬ 
sistent with the presence of a first-order phase bound¬ 
ary between the amorphous phases of ice, as noted in 
Ref. 100,101 . Furthermore, this finding complements the 
distinctive features of HDA and LDA observed via the 
goo{r) above. 


D. Autocorrelation of the LSI Order Parameter in 
Ambient Liquid Water 


All of the analyses and discussion provided above in 
Section III C were only made possible with the knowledge 
of the liquid water structures at the IPES, obtained by 
removing the thermal motion effects from ambient liq¬ 
uid water.At ambient temperatures, these thermal mo¬ 
tions smear the distinctions made based on in the 
IPES and effectively result in the unimodal P(J) dis¬ 
cussed in Section III A (see Figure 2). Hence, these ther¬ 
mal motions introduce decay in the fluctuations in the 
local structural environment surrounding a given water 
molecule. 

To investigate the temporal evolution of these thermal 
fluctuations in the local structural environment, we have 
computed the autocorrelation function, C/(i), of the LSI 
order parameter, defined as: 


(61(0)61®) 
n ’ ( 61(0)61(0 )) 


( 2 ) 


in which 61(0) ( 6I(t )) is defined for a given water 
molecule as the difference between I given by Equa¬ 
tion (1) at time t = 0 (t = t) and the ensemble average 
I over all molecules and MD snapshots, i.e., 61 = I — I. 
As shown in Figure 6, the LSI autocorrelation function 
decays to zero within approximately 4 ps for the case 
of PBE0+TS-vdW(SC) liquid water computed at 330 K. 
Here we note that this ~ 4 ps persistence time corre¬ 
sponding to different local molecular environments is sim¬ 
ilar to the time scales associated with the density-density 




























IV. CONCLUSIONS AND FUTURE OUTLOOK 
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FIG. 6. Autocorrelation function of the local structure index 
order parameter, Ci(t ), computed from DFT-based AIMD 
(PBEO+TS-vdW(SC)) and classical force field (TIP4P/Ice) 
based simulations of liquid water at temperatures in the range 
of 300-330 K. 



fluctuations observed in large-scale classical simulations 
of ambient liquid water 38 . For comparison, the C/(f) 
obtained from TIP4P/Ice liquid water at 330 K shows 
almost identical behavior with the AIMD results. Upon 
reduction of the simulation temperature to 300 K, the 
autocorrelation functions for both PBE0+TS-vdW(SC) 
and TIP4P/Ice still remain extremely similar with cor¬ 
relation times increasing to ~ 6 ps. 

From Figure 6, one can also notice that the decay rate 
of Ci(t) below 0.5 ps occurs significantly faster and is es¬ 
sentially independent of the simulation temperature—an 
effect which is primarily due to the librational motion in 
liquid water. Beyond 0.5 ps, the decay rate of Cj(t) is 
much slower and is influenced by the collective thermal 
motions of the molecules—motions which are fairly sensi¬ 
tive to the simulation temperature. These findings show 
that naturally occurring thermal fluctuations in ambient 
liquid water originate local LSI fluctuations which allow 
water molecules to sample different local structural envi¬ 
ronments. Hence the presence of low- and high-density¬ 
like sites are rather transient, with a given molecule sam¬ 
pling the continuous distribution of local environments 
existing between these two distinct states over the course 
of approximately ~ 4 ps. 

In addition, we note in passing that the sluggish dy¬ 
namics observed in liquid water simulated at 300 K uti¬ 
lizing the PBE XC functional are directly reflected in the 
persistence time between low- and high-density-like sites. 
The Ciit) computed from PBE (300 K) liquid water de¬ 
cays at a much slower rate and vanishes at ~ 14 ps, which 
is approximately 3 times longer than that computed from 
PBE0+TS-vdW(SC) (330 K) liquid water. 


In this work, we have performed a detailed computa¬ 
tional study of the fluctuations in the local structural 
environments found in ambient liquid water using the 
local structure index (LSI) order parameter and a sys¬ 
tematic hierarchy of DFT XC functionals. For each of 
the functionals considered herein, we found probability 
distributions of the LSI order parameter, P(/), that were 
all unimodal in shape and rapidly decaying from max¬ 
ima located in the low-LSI range, indicating that a ma¬ 
jority of the water molecules at ambient conditions are 
situated in locally disordered high-density-like environ¬ 
ments. We observed large variations among the relative 
populations of the continuous distribution of local struc¬ 
tural environments existing in ambient liquid water as a 
function of the underlying XC potential —a trend that 
appears to be remarkably similar to the effects observed 
by changing the simulation temperature. At the PBE 
level of theory, there exists a much higher percentage of 
molecules with high-LSI values, ie., molecules which are 
characterized by low-density-like and more ordered local 
environments, and this finding provides additional mi¬ 
croscopic structural insight into the overstructuring and 
sluggish nature of PBE liquid water simulated at ambi¬ 
ent conditions. Including the effects of exact exchange 
and vdW/dispersion interactions in the underlying XC 
potential resulted in a significant increase in the relative 
population of molecular sites with low-LSI values, which 
are more disordered and high-density like. This increased 
population of high-density-like molecules originates from 
an energetic weakening of the hydrogen bonds in conjunc¬ 
tion with an entropic stabilization of the resultant disor¬ 
dered water configurations, as facilitated by the inclusion 
of exact exchange and non-local vdW/dispersion forces. 
As a result, ambient liquid water simulated with the 
vdW-inclusive hybrid functional, PBE0+TS-vdW(SC), 
is significantly less structured and more diffusive than 
PBE water with a substantial improvement in the agree¬ 
ment with experiment for a number of structural and 
electrostatic properties 1 . 

Focusing on liquid water generated at the PBE0+TS- 
vdW(SC) level of theory, we found the P(7) to be uni¬ 
modal at ambient temperatures, with removal of ther¬ 
mal excitation leading to a bimodal P(J) distribution at 
the corresponding IPES. This behavior is reflected in the 
density and spatial distribution of the low-density-like 
molecules. In the IPES, the population of low-density¬ 
like molecules grows considerably and is accompanied 
by the formation of connected clusters which can in¬ 
clude as many as 34 molecules in a (H20)i28 simula¬ 
tion cell. These clusters do not span the length of the 
cell as found at liquid-liquid coexistence 23 , but are suf¬ 
ficiently extended to allow for meaningful ring statistics 
analysis. The presence of these low-density-like clusters 
amid the larger population of high-density-like molecules 
in the IPES of ambient water is noteworthy, as these 
clusters should act like seeds for the growth of larger low- 
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density clusters upon cooling. This behavior has indeed 
been observed in MD simulations with empirical force 
fields upon cooling 33-35 and is entirely consistent with 
the experimental observation of metastable /mi-sized liq¬ 
uid droplets with low-density character 22 . It is the pres¬ 
ence of low-density and high-density fluctuations grow¬ 
ing in size upon cooling that leads, under appropriate 
thermodynamic and kinetic conditions, to the formation 
of the LDA and/or HD A glasses and their liquid coun¬ 
terparts, LDL and/or HDL. Our finding that the sig¬ 
natures of the intermediate range order (IRO) of LDA 
and HDA are both present in the IPES of ambient wa¬ 
ter is therefore quite remarkable. It is also interesting 
that these signatures are associated to ring statistics dis¬ 
tributions that correspond to distinct topologies of the 
underlying hydrogen bond network. In this regard, a 
change of topology from low-density-like to high-density¬ 
like involves substantial breaking and reforming of hydro¬ 
gen bonds accompanied by an associated energy barrier, 
which is consistent with the presence of a first-order tran¬ 
sition boundary between LDA and HDA 100,101 . 

Analysis of the autocorrelation of the LSI order pa¬ 
rameter at ambient conditions shows that low- and high- 
density fluctuations are transient and decay with a life¬ 
time of ~ 4 ps at the PBE0+TS-vdW(SC) level of theory. 
In other words, the naturally occurring thermal fluctua¬ 
tions allow a given water molecule to sample the continu¬ 
ous distribution of local structural environments between 
the two distinct high- and low-density states every ~ 4 
ps on average. In this respect, it should also be noted 
that alternative two-state hydrogen bond models have 
often been invoked to explain the features observed in 
experiments probing the local structure of ambient wa- 

j- er 29, 30,104-106 

As a final remark, we should mention that the quan¬ 
titative details of the picture that we have presented in 
this work may change by further refining the ab initio 
description of water with more accurate XC functionals 


that reduce the self-interaction error, as could possibly 
be achieved by fine-tuning the exchange correction 107 , 
and/or include a better description of vdW/dispersion 
interactions via the inclusion of beyond-pairwise interac¬ 
tions, as provided, for example, by the recently proposed 
many-body dispersion (MBD) scheme 108-111 . Further re¬ 
finements are expected to result from a proper treatment 
of nuclear quantum effects via explicit treatment of the 
quantum mechanical nature of the nuclei by the Feyn¬ 
man discretized path-integral approach. However, all of 
these refinements are unlikely to change the qualitative 
picture exposed herein. 
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